Practical 2 - hospital admissions for respiratory conditions in Greater Glasgow and Clyde (2007-2011)

In this practical you will use R (R Core Team 2019) as well as INLA to carry out a spatial and a spatio-temporal small area disease risk analysis.

In particular, you are going to model yearly hospital admissions for respiratory conditions (ICD-10 codes: J00-J99) in the Greater Glasgow and Clyde health board, for the period 2007 - 2011.

Scotland is divided into 14 health boards, and here we focus on the Greater Glasgow and Clyde health board, which contains the city of Glasgow and has a population of around 1.2 million people, during the period 2007 to 2011. This health board is split into N = 271 intermediate geographies (IG), which are a key geography for the distribution of small-area statistics in Scotland and contain populations of between 2,468 and 9,517 people.

The data used in this practical are freely available from the R package CARBayesST (Lee, Rushworth, and Napier 2018). A subset of these data is in the file RESP_DATA.csv.

Before starting the practical

• Load needed libraries:

library(dplyr)        # A package for data manipulation
library(sf)           # Simple feature for R
library(spdep)        # Functions and tests for evaluating spatial patterns 
# and autocorrelation
library(tidyr)

library(INLA)         # Integrated Nested Laplace Approximation package
library(ggplot2)      # A package that implements the grammar of graphics, which is a term used to
# break up graphs into semantic components, such as geometries and layers.
library(viridis)      # A package providing color palettes 
library(patchwork)

# For tables in RMarkdown
library(knitr)
library(kableExtra)

1. Data

1. Import the.csv file with the data and call the data.frame object as RESP_DATA.

RESP_DATA <- read.csv("RESP_DATA.csv", header=TRUE)

Here The first column labelled IG is the set of unique identifiers for each IG, year is the year of hospitalization, while observed and expected are respectively the observed and expected numbers of hospital admissions due to respiratory diseases.

2. Compute the total number of cases of hospital admissions per year, and format the output in a table.

kable(RESP_DATA %>%
        group_by(year) %>%
         summarise(observed = sum(observed), expected=sum(expected)), booktabs = T, caption = "Hospital admissions by year") %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "center")

Hospital admissions by year

year	observed	expected
2007	20410	23727.31
2008	21960	24955.70
2009	21166	24206.05
2010	21234	26007.28
2011	22548	26234.11

Note: to know more about this function, you can type in the console ?knitr::kable and visit the page https://haozhu233.github.io/kableExtra/awesome_table_in_pdf.pdf

3. Import the shape file GGHB using the function st_read from sf package and call the object as GGHB.

GGHB <- st_read("GGHB.shp")

## Reading layer `GGHB' from data source 
##   `/Users/marta/Dropbox/Books/INLABook/ShortCourse/VIBASS/Practicals/Practical2/GGHB.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 271 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 218820 ymin: 644792 xmax: 272283 ymax: 689565
## CRS:           NA

4. Then, plot the spatial object GGHB using ggplot2 package.

ggplot() + 
      geom_sf(data = GGHB, color = "blue", fill = "white") + 
      coord_sf() +    #axis limits and CRS
      theme_bw() +    # dark-on-light theme
      theme(axis.title = element_text(size = 14),
            axis.text = element_text(size = 12))

The figure shows that the river Clyde partitions the study region into a northern and a southern sub-region, and no areal units on opposite banks of the river border each other. We will account for this issue when creating neighbors for the areal units.

2. Spatial model

In order to use the same data structure for both the space-only model and later the space-time model, a new set of data is formed by aggregating both the observed and expected counts over time. A Poisson-log linear model is then fitted, assuming a BYM2 model for the random effects. Let each areal unit \(i\) be indexed by the integers \(1, 2,...,N\). \[ \begin{eqnarray} O_{i}|\lambda_{i} & \sim & \text{Poisson}(\lambda_{i}E_{i} ) \\ \log(\lambda_{i}) & = & \alpha + \theta_{i} + \phi_{i} \\ \theta_{i} & \sim & N(0, \sigma_{\theta}^2)\\ {\bf \phi} & \sim & \hbox{ICAR}({\bf W}, \sigma_{\phi}^2)\\ \alpha & \sim & \text{Uniform}(-\infty, +\infty) \\ 1/\sigma_{\theta}^2 & \sim & \text{Gamma}(1,0.001) \\ 1/\sigma_{\phi}^2 & \sim & \text{Gamma}(0.5,0.005) \\ \end{eqnarray} \]

5. Define the neighbors and create the weights list. Due to the river, some areas are not connected. To avoid this artifact you can use the snap argument (boundary points less than snap distance apart are considered to indicate contiguity; see https://cran.r-project.org/web/packages/spdep/vignettes/nb.pdf)

GGHB_nb <- poly2nb(GGHB, snap=1000, queen=TRUE)
summary(GGHB_nb)

## Neighbour list object:
## Number of regions: 271 
## Number of nonzero links: 3414 
## Percentage nonzero weights: 4.64863 
## Average number of links: 12.59779 
## Link number distribution:
## 
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 23 25 27 33 
##  1  4  6  7  6 11 16 12 17  8 18 15 27 21 23 18 19 17  7  6  4  5  1  1  1 
## 1 least connected region:
## 202 with 1 link
## 1 most connected region:
## 234 with 33 links

Convert the list of neighbors to inla format using the function nb2WB().

nb2INLA("GGHB.graph",GGHB_nb)
GGHB.adj <- paste(getwd(),"/GGHB.graph",sep="")

6. Aggregate observed and expected cases over geographical areas

RESP_DATA %>% group_by(SP_ID) %>% 
              summarize(observed = sum(observed), 
                        expected = sum(expected)) %>% 
              dplyr::rename(O = observed, E = expected) -> RESP_DATAagg

7. Compute the standardized morbidity ratios (SMRs)

RESP_DATAagg %>% mutate(SMR = O/E) -> RESP_DATAagg

8. Produce a spatial map of the aggregated SMRs using ggplot2 package. For the map use the following breakpoints [min,0.4], (0.4-0.6], (0.6-0.8], (0.8,1], (1,1.2], (1.2-1.4], (1.4-1.6], (1.6-max].

Remember that, before to produce the map, you need to join the sf object and the data frameRESP_DATAagg. To do so you can use the function left_join from the library dplyr (see previous practicals).

RESP_DATAagg$SMRcat <- cut(RESP_DATAagg$SMR, 
                      breaks=c(min(RESP_DATAagg$SMR), 
                               0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 
                               max(RESP_DATAagg$SMR)), include.lowest = T)

map_SMR <- left_join(GGHB, RESP_DATAagg, by = c("SP_ID" = "SP_ID"))

and plot:

ggplot() + geom_sf(data = map_SMR, col = NA) + aes(fill = SMRcat) +
  theme_bw() + scale_fill_viridis_d() + 
  guides(fill=guide_legend(title="SMR")) 

Map of the average SMRs over the period 2007-2011

9. Fit the hierarchical Poisson log-linear model in INLA

ID<- seq(1,271)
formula_BYM2 <- O ~ f(ID, model="bym2", graph=GGHB.adj,
                            hyper=list(prec = list(
        prior = "pc.prec",
        param = c(0.5 / 0.31, 0.01)),
        phi = list(
        prior = "pc",
        param = c(0.5, 2 / 3))))    
sBYM.model <- inla(formula=formula_BYM2, family="poisson", data=RESP_DATAagg, E=E, control.compute=list(dic=TRUE, waic=TRUE))

10. Obtain the posterior summary statistics (mean and posterior probability that the residual is above 1 - (or log-residual is above 0)) of the parameters of interest

#Relative risks
RR_sBYM<-c()
for(i in 1:271){
  RR_sBYM[i] <- inla.emarginal(function(x) exp(x), 
        sBYM.model$marginals.random$ID[[i]])
}

#Posterior probabilities
RR_sBYM_marg <- sBYM.model$marginals.random$ID[1:271]
PP_sBYM <- lapply(RR_sBYM_marg, function(x) {1-inla.pmarginal(0,x)})

11. Obtain the posterior estimates from the spatial model to be plotted, that is (i) the area level posterior mean of the residual RRs and (ii) the posterior probability (PP) that the residual RRs > 1.

resRR_PP <- data.frame(resRR=RR_sBYM, 
                       PP=unlist(PP_sBYM),
                      SP_ID=RESP_DATAagg[,1])

16. Using ggplot2 package, produce a map of the posterior mean of the residual RRs and the posterior probabilities that the residual RRs are > 1

• For the map of the posterior mean of the residual RRs, use the following breakpoints [min,0.4], (0.4-0.6], (0.6-0.8], (0.8,1], (1,1.2], (1.2-1.4], (1.4-1.6], (1.6-max].

resRR_PP$resRRcat <- cut(resRR_PP$resRR, breaks=c(min(resRR_PP$resRR), 
                  0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 
                  max(resRR_PP$resRR)),include.lowest = T)

• For the map of the probabilities PP that the residual RRs is > 1, use the following breakpoints [0-0.2], (0.2-0.8], (0.8-1].

# breakpoints
resRR_PP$PPcat <- cut(resRR_PP$PP, c(0, 0.2, 0.8, 1.00), include.lowest = TRUE)

• Remember to join sf object and data frame with the posterior estimates

map_RR_PP <- left_join(GGHB, resRR_PP, by = c("SP_ID" = "SP_ID"))

• Produce the maps of the posterior mean of the residual RRs and the posterior probabilities PP using ggplot2 package.

ggplot() + geom_sf(data = map_RR_PP) + aes(fill = resRRcat) +
  theme_bw() + scale_fill_brewer(palette = "PuOr") + 
  guides(fill=guide_legend(title="RR")) + ggtitle("RR Spatial model") + 
  theme(text = element_text(size=15), 
                  axis.text.x = element_blank(), 
                  axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold")
                  )-> p1

ggplot() + geom_sf(data = map_RR_PP) + aes(fill = PPcat) +
  theme_bw() +
  scale_fill_viridis(
    option = "plasma", name="PP",
    discrete = T,
    direction = -1,
    guide = guide_legend(
      title.position = 'top',
      reverse = T
    )) +  ggtitle("PP Spatial model") + theme(text = element_text(size=15), 
                  axis.text.x = element_blank(), 
                  axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold")
                  ) -> p2

p1|p2

3. Estimate the spatial fraction

As the BYM2 has the structured (CAR) and unstructured (iid) components it might be useful to get some ideas about the strength of the spatially structured components as this would indicate the level of clustering in the data. To do so we ca simply obtain the posterior summary of the phi hyperparameter

sBYM.model$summary.hyperpar

##                       mean        sd 0.025quant  0.5quant 0.975quant mode
## Precision for ID 7.6801757 0.8105148  6.1734660 7.6512274   9.362810   NA
## Phi for ID       0.6716711 0.1257730  0.4069953 0.6802351   0.887077   NA

This tells us that about 2/3 of the spatial variability is explained by the spatially structured component - which makes sense if we look at the map of the RR, which show a degree of spatial clustering.

4. Spatio-temporal model (no interaction)

Now, we extend the above analysis to a separable space-time model without interactions. For the temporal component, we use the specification introduces in the lecture with a temporal unstructured random effect and a structured one (RW1 prior).

Let each areal unit \(i\) be indexed by the integers \(1, 2,...,N\). As in the spatial case, we use a Poisson distribution to model the number of hospital admission \(O_it\), in area \(i\) at time \(t\). The mathematical specification of the model includes now an additional temporal dependence term, which can be modeled using a non-stationary random walk prior: ${i} ({i-1}, ^2_{}).
The model implement in this practical assumes no space-time interaction and a spatial convolution with random walk in time:

\[ \begin{eqnarray} O_{i}|\lambda_{i} & \sim & \text{Poisson}(\lambda_{i}E_{i} ) \\ \log(\lambda_{i}) & = & \alpha + \theta_{i} + \phi_{i} + \xi_t + \gamma_t \\ \theta_{i} & \sim & N(0, \sigma_{\theta}^2)\\ {\bf \phi} & \sim & \hbox{ICAR}({\bf W}, \sigma_{\phi}^2)\\ \gamma_t &\sim & N(0, \sigma_{\gamma}^2) \\ \xi_t & \sim & RW1(\sigma_{\xi}^2) \\ \alpha & \sim & \text{Uniform}(-\infty, +\infty) \\ 1/\sigma_{\theta}^2 & \sim & \text{Gamma}(1,0.001) \\ 1/\sigma_{\phi}^2 & \sim & \text{Gamma}(0.5,0.005) \\ 1/\sigma_{\gamma}^2 & \sim & \text{Gamma}(1,0.001) \\ 1/\sigma_{\xi}^2 & \sim & \text{Gamma}(0.5,0.005) \\ \end{eqnarray} \]

17. First prepare the data, joining in the shapefile to make sure that the order is the same and then create an ID for time and one for space.

#Join the data with the shapefile so the order of the shapefile is maintained.  
RESP_DATA_ST <- left_join(GGHB, RESP_DATA, by="SP_ID")
#Rename the columns of Observed and Expected as we did before
RESP_DATA_ST <- RESP_DATA_ST  %>% dplyr::rename(O = observed, E = expected)
#Create the ID for year (time)
RESP_DATA_ST$ID.time <- RESP_DATA_ST$year - 2006
#Create the ID for space
RESP_DATA_ST$ID.space <- rep(seq(1,271),each=5)

Run the model in INLA

formula_ST_noint <- O ~ f(ID.space, model="bym2", graph=GGHB.adj,
                            hyper=list(prec = list(
                            prior = "pc.prec",
                            param = c(0.5 / 0.31, 0.01)),
                            phi = list(
                            prior = "pc",
                            param = c(0.5, 2 / 3)))) + f(ID.time,model="rw1", hyper=list(prec = list(
                            prior = "pc.prec",
                            param = c(0.5 / 0.31, 0.01))))
                            
stBYM.model <- inla(formula=formula_ST_noint, family="poisson", data=RESP_DATA_ST, E=E, control.compute=list(dic=TRUE, waic=TRUE))

18. Create the posterior mean for the spatial (as we did in point 10) and temporal effects

#Spatial Relative risks
RR_stBYM<-c()
for(i in 1:271){
  RR_stBYM[i] <- inla.emarginal(function(x) exp(x), 
        stBYM.model$marginals.random$ID.space[[i]])
}
#Posterior probabilities (for spatial RR)
RR_stBYM_marg <- stBYM.model$marginals.random$ID.space[1:271]
PP_stBYM <- lapply(RR_stBYM_marg, function(x) {1-inla.pmarginal(0,x)})

#Temporal Relative risks and CI95
RR_stRW_RR<-c()
RR_stRW_lo<-c()
RR_stRW_hi<-c()

for(i in 1:5){
  #Posterior mean
  RR_stRW_RR[i] <- inla.emarginal(function(x) exp(x), 
        stBYM.model$marginals.random$ID.time[[i]])
  #2.5% quantile 
  RR_stRW_lo[i] <- inla.qmarginal(0.025,inla.tmarginal(function(x) exp(x), stBYM.model$marginals.random$ID.time[[i]]))
  #97.5% quantile 
  RR_stRW_hi[i] <- inla.qmarginal(0.975, inla.tmarginal(function(x) exp(x), stBYM.model$marginals.random$ID.time[[i]]))
}

RR_stRW <- data.frame(RR=RR_stRW_RR,low=RR_stRW_lo,high=RR_stRW_hi)

19. Plot the temporal residual RRs (RR_stWR)

ggplot(RR_stRW, aes(seq(2007,2011), RR)) + geom_line() + ggtitle("ST model No Int") + geom_ribbon(aes(ymin=low,ymax=high), alpha=0.2) + labs(x="year")-> Temp1
Temp1

20. Map the spatial residual RRs (RR_stBYM) with ggplot2 package using the following breakpoints [min,0.4], (0.4-0.6], (0.6-0.8], (0.8,1], (1,1.2], (1.2-1.4], (1.4-1.6], (1.6-max]. Compare this map against the map of the residual RR obtained from the spatial model.

resRR_PP_st <- data.frame(resRR=RR_stBYM, 
                       PP=unlist(PP_stBYM),
                      SP_ID=RESP_DATAagg[,1])
# breakpoints
resRR_PP_st$resRRcat <- cut(resRR_PP_st$resRR, breaks=c(min(resRR_PP_st$resRR), 
                  0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 
                  max(resRR_PP_st$resRR)),include.lowest = T)

resRR_PP_st$PPcat <- cut(resRR_PP_st$PP, c(0, 0.2, 0.8, 1.00), include.lowest = TRUE)

map_RR_ST <- left_join(GGHB, resRR_PP_st, by = c("SP_ID" = "SP_ID"))

ggplot() + geom_sf(data = map_RR_ST) + aes(fill = resRRcat) +
  theme_bw() + scale_fill_brewer(palette = "PuOr") + 
  guides(fill=guide_legend(title="RR")) +  ggtitle("RR ST model") +
  theme(text = element_text(size=15), 
        axis.text.x = element_blank(), 
        axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold")
        ) -> p3

ggplot() + geom_sf(data = map_RR_ST) + aes(fill = PPcat) +
  theme_bw() +
  scale_fill_viridis(
    option = "plasma",
    name = "PP ST model",
    discrete = T,
    direction = -1,
    guide = guide_legend(
      title.position = 'top',
      reverse = T
    )) +  ggtitle("PP ST model") + theme(text = element_text(size=15), 
        axis.text.x = element_blank(), 
        axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold")
        )-> p4

(p1|p2) / (p3|p4)

Spatio-temporal model: Map of the residual RRs and posterior probabilities

5. Spatio-temporal model (type I interaction)

Now, we extend the above analysis to a separable space-time model with type I interaction:

\[ \begin{eqnarray} O_{i}|\lambda_{i} & \sim & \text{Poisson}(\lambda_{i}E_{i} ) \\ \log(\lambda_{i}) & = & \alpha + \theta_{i} + \phi_{i} + \xi_t + \gamma_t + \zeta_{it} \\ \theta_{i} & \sim & N(0, \sigma_{\theta}^2)\\ {\bf \phi} & \sim & \hbox{ICAR}({\bf W}, \sigma_{\phi}^2)\\ \gamma_t &\sim & N(0, \sigma_{\gamma}^2) \\ \xi_t & \sim & RW1(\sigma_{\xi}^2) \\ \zeta_{it} & \sim & N(0, \sigma_{\zeta}^2) \\ \alpha & \sim & \text{Uniform}(-\infty, +\infty) \\ 1/\sigma_{\theta}^2 & \sim & \text{Gamma}(1,0.001) \\ 1/\sigma_{\phi}^2 & \sim & \text{Gamma}(0.5,0.005) \\ 1/\sigma_{\gamma}^2 & \sim & \text{Gamma}(1,0.001) \\ 1/\sigma_{\xi}^2 & \sim & \text{Gamma}(0.5,0.005) \\ 1/\sigma_{\zeta}^2 & \sim & \text{Gamma}(1,0.001) \end{eqnarray} \]

21. Specify the formula and run the model in INLA. Remember that you need to create an index which goes from 1 to the length of the dataset (i.e. the space x time)

RESP_DATA_ST$ID.space.time <- seq(1,dim(RESP_DATA_ST)[1])
formula_ST_intI <- O ~ f(ID.space, model="bym2", graph=GGHB.adj,
                            hyper=list(prec = list(
                            prior = "pc.prec",
                            param = c(0.5 / 0.31, 0.01)),
                            phi = list(
                            prior = "pc",
                            param = c(0.5, 2 / 3)))) + 
                      f(ID.time,model="rw1", hyper=list(prec = list(
                            prior = "pc.prec",
                            param = c(0.5 / 0.31, 0.01))))+
                      f(ID.space.time,model="iid", hyper=list(prec = list(
                            prior = "pc.prec",
                            param = c(0.5 / 0.31, 0.01))))
                    
                            
stIntI.BYM.model <- inla(formula=formula_ST_intI, family="poisson", data=RESP_DATA_ST, E=E, control.compute=list(dic=TRUE, waic=TRUE))

22. Create the posterior mean for the spatial and temporal effects and compare with the ST model results without interaction

#Spatial Relative risks
RR_stIntI.BYM<-c()
for(i in 1:271){
  RR_stIntI.BYM[i] <- inla.emarginal(function(x) exp(x), 
        stIntI.BYM.model$marginals.random$ID.space[[i]])
}
#Posterior probabilities (for spatial RR)
RR_stIntI.BYM_marg <- stIntI.BYM.model$marginals.random$ID.space[1:271]
PP_stIntI.BYM <- lapply(RR_stIntI.BYM_marg, function(x) {1-inla.pmarginal(0,x)})

#Temporal Relative risks and CI95
RR_stIntI.RW_RR<-c()
RR_stIntI.RW_lo<-c()
RR_stIntI.RW_hi<-c()

for(i in 1:5){
  #Posterior mean
  RR_stIntI.RW_RR[i] <- inla.emarginal(function(x) exp(x), 
        stIntI.BYM.model$marginals.random$ID.time[[i]])
  #2.5% quantile 
  RR_stIntI.RW_lo[i] <- inla.qmarginal(0.025,inla.tmarginal(function(x) exp(x), stIntI.BYM.model$marginals.random$ID.time[[i]]))
  #97.5% quantile 
  RR_stIntI.RW_hi[i] <- inla.qmarginal(0.975, inla.tmarginal(function(x) exp(x), stIntI.BYM.model$marginals.random$ID.time[[i]]))
}

RR_stIntI.RW<- data.frame(RR=RR_stIntI.RW_RR,low=RR_stIntI.RW_lo,high=RR_stIntI.RW_hi)

23. Plot the temporal residual RRs (RR_stWR)

ggplot(RR_stIntI.RW, aes(seq(2007,2011), RR)) + geom_line() + ggtitle("ST model Int I") + geom_ribbon(aes(ymin=low,ymax=high), alpha=0.2) + labs(x="year")->Temp2
Temp1 | Temp2

24. Map the spatial residual RRs (RR_stIntI.BYM) with ggplot2 package using the following breakpoints [min,0.4], (0.4-0.6], (0.6-0.8], (0.8,1], (1,1.2], (1.2-1.4], (1.4-1.6], (1.6-max]. Compare this map against the map of the residual RR obtained from the spatial temporal model with no interaction.

resRR_PP_stIntI <- data.frame(resRR=RR_stIntI.BYM, 
                       PP=unlist(PP_stIntI.BYM),
                      SP_ID=RESP_DATAagg[,1])
# breakpoints
resRR_PP_stIntI$resRRcat <- cut(resRR_PP_stIntI$resRR, breaks=c(min(resRR_PP_stIntI$resRR), 
                  0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, 
                  max(resRR_PP_stIntI$resRR)),include.lowest = T)

resRR_PP_stIntI$PPcat <- cut(resRR_PP_stIntI$PP, c(0, 0.2, 0.8, 1.00), include.lowest = TRUE)

map_RR_ST.IntI <- left_join(GGHB, resRR_PP_stIntI, by = c("SP_ID" = "SP_ID"))

ggplot() + geom_sf(data = map_RR_ST.IntI) + aes(fill = resRRcat) +
  theme_bw() + scale_fill_brewer(palette = "PuOr") + 
  guides(fill=guide_legend(title="RR")) +  ggtitle("RR ST model Int I") +
  theme(text = element_text(size=15), 
        axis.text.x = element_blank(), 
        axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold")
        ) -> p5

ggplot() + geom_sf(data = map_RR_ST.IntI) + aes(fill = PPcat) +
  theme_bw() +
  scale_fill_viridis(
    option = "plasma",
    name = "PP ST model Int I",
    discrete = T,
    direction = -1,
    guide = guide_legend(
      title.position = 'top',
      reverse = T
    )) +  ggtitle("PP ST model Int I") + theme(text = element_text(size=15), 
        axis.text.x = element_blank(), 
        axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold")
        )-> p6

(p1|p2) / (p3|p4) / (p5|p6)

Spatio-temporal model: Map of the residual RRs and posterior probabilities

We basically see that across the different models there is no difference in the spatial residuals. Let’s now look at the ST interaction.

25. Plot the spacetime interaction.

RESP_DATA_ST$intI<-stIntI.BYM.model$summary.random$ID.space.time$mean
RESP_DATA_ST$intI_cat <- cut(RESP_DATA_ST$intI,  breaks=c(-1,-0.05, 
                  -0.01, 0.01, 0.05, 1),include.lowest = T)
ggplot() +
  geom_sf(data = RESP_DATA_ST, aes(fill = intI_cat))+ theme_bw() +  scale_fill_brewer(palette = "PuOr") + 
            guides(fill=guide_legend(title=NULL)) + 
            theme(text = element_text(size=20), 
                  axis.text.x = element_blank(), 
                  axis.text.y = element_blank()) +
facet_wrap(~ year, ncol = 3, labeller=labeller(ID.year=c("1"="2007","2"="2008","3"="2009","4"="2010","5"="2011"))) +
labs("")

We can see that there is not clear pattern in the interactions.

• Get a table of the hyperparameters. How do you interpret this table?

dat.hyper2 <- 
  round(
  data.frame(median = stIntI.BYM.model$summary.hyperpar[,4],
    LL = stIntI.BYM.model$summary.hyperpar[,3], 
    UL = stIntI.BYM.model$summary.hyperpar[,5]),
  digits = 3)

row.names(dat.hyper2) <- 
  rownames(stIntI.BYM.model$summary.hyperpar)

knitr::kable(dat.hyper2, caption = "Posterior median and 95% CrI of hyperparameters.") %>%  kable_styling(bootstrap_options = "striped", full_width = F, position = "center")

Posterior median and 95% CrI of hyperparameters.

	median	LL	UL
Precision for ID.space	7.862	6.277	9.628
Phi for ID.space	0.680	0.394	0.911
Precision for ID.time	362.430	72.249	2122.259
Precision for ID.space.time	117.504	95.447	146.324

• Compare the WAIC. What do you observe?

dat.WAIC <- data.frame(model = c("Spatial", "SpatTemp no int", "SpatTemp typeI"), 
                       WAIC = round(c(sBYM.model$waic$waic, stBYM.model$waic$waic, stIntI.BYM.model$waic$waic))
)

row.names(dat.WAIC) <- NULL

knitr::kable(dat.WAIC, caption = "WAIC of the fifferent models") %>%  kable_styling(bootstrap_options = "striped", full_width = F, position = "center")

WAIC of the fifferent models

model	WAIC
Spatial	2558
SpatTemp no int	10820
SpatTemp typeI	10312

Lee, D., A. Rushworth, and G. Napier. 2018. “Spatio-Temporal Areal Unit Modeling in r with Conditional Autoregressive Priors Using the CARBayesST Package.” Journal of Statistical Software, Articles 84 (9).

R Core Team. 2019. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.